Dual Fermion Approach to Susceptibility of Correlated Lattice Fermions 
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In this paper, we show how the two-particle Green lunction (2PGF) can be obtained within 
the Iramework ol the Dual Fermion approach. This lacilitates the calculation ol the susceptibility 
in strongly correlated systems where long-ranged non-local correlations cannot be neglected. We 
formulate the Bethe-Salpeter equations for the full vertex in the particle-particle and particle-hole 
channels and introduce an approximation for practical calculations. The scheme is applied to the 
two-dimensional Hubbard model at half filling. The spin-spin susceptibility is found to strongly 
increase for the wavevector q = (n,iv), indicating the antiferromagnetic instability. We find a 
suppression of the critical temperature compared to the mean-field result due to the incorporation 
of the non-local spin-fluctuations. 
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I. INTRODUCTION 

Strongly correlated electron systems exhibit some of 
the most intriguing features known to condensed mat- 
ter physics, including high-temperature superconductiv- 
ity, heavy-fermion behavior, different kinds of electronic 
phase transitions, etcJ*&&^&. In most of the cases, a 
conventional band theory 6 '- is not sufficient to describe 
these properties. Essentially a many-body treatment is 
necessary, correlation effects being too strong to be de- 
scribed by the standard perturbation theory. Due to the 
complexity of this problem the physics of strongly cor- 
related systems has proven to be theoretically challeng- 
ing. One of the successful routes to the description of 
strongly correlated systems is the Dynamical Mean Field 
Theory (DMFT) 8 ' 9 . It is commonly accepted now that 
this approach typically catches the most essential cor- 
relation effects, e.g., the physics of the Mott-Hubbard 
transition^ 9 -. The method was implemented success- 
fully into realistic electronic structure calculations 1 ^!!, 
which is now a standard tool in the microscopic theory of 
strongly correlated systems^ 2 -. However, there are many 
phenomena for which non-local correlations are impor- 
tant and often the relevant correlations are even long- 
ranged. Among them are the Luttinger-liquid formation 
in low-dimensional systems&i^, non-Fermi-liquid behav- 
ior due to van-Hove singularities in two dimension o 14 ' 15 , 
the physics near quantum critical points 16 , or d-wave 
pairing in high-T c superconductors^. 

The requirement to incorporate spatial correlations 
of strongly correlated fermions into present theories 
has triggered several efforts to go beyond DMFT. In 
the Dynamical Vertex Approximation 1 ^ and similar 
approache s 18 ' 19 , a diagrammatic expansion around the 
DMFT solution is performed. In Ref. the authors 
introduce a scheme based on the assumption of the local- 
ity of the fully irreducible vertex of the lattice problem. 
In the framework of their approximation this vertex is 



extracted from the Anderson impurity model. The re- 
ducible vertex can be obtained from the parquet equa- 
tions which is subsequently used to calculate the non- 
local self-energy. The Green function is updated and the 
parquet equations are solved again until self-consistency 
is reached. In the work cited above a simplified version 
of this algorithm was implemented, where instead of the 
parquet equations the Bethe-Salpeter equation (BSE) in 
a particular channel was solved without performing a self- 
consistent calculation. Hence, this approach does not go 
beyond the usual DMFT approximation for calculating 
the two-particle Green functions (2PGFs). Practically, to 
go well beyond DMFT within this approach, the parquet 
equations, which form a complicated system of coupled 
integral equations need to be solved repeatedly, which 
requires a sizeable numerical effort. 

A principally new scheme with a fully renormalizcd 
expansion called Dual Fermion Approach has been pro- 
posed recently 20 . It is based on the introduction of new 
variables in the path integral representation. This ap- 
proach yields very satisfactory results already for the 
lowest-order corrections, while the schemes proposed 
in Refsj 17 ' 18 ' 19 operate with infinite diagrammatic se- 
ries and require the solution of complicated integral 
equations. A scheme similar to the Dual Fermion ap- 
proach has been discussed earlier in terms of Hubbard 
operators^, but has not been used for practical calcula- 
tions. 

In many cases, it is desirable to calculate the 2PGF 
which provides insight into the inter-particle interactions 
and two-particle excitations of the system and corre- 
sponding instabilities. In the framework of DMFT the 
2PGFs are usually calculated under the assumption that 
the irreducible vertex part in a given channel is local and 
it is taken to be equal to the corresponding irreducible 
vertex of the impurity problem. This assumption is em- 
pirical and does not hold in the general case 2 ^. 

In this paper we describe how to perform calculations 
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of the 2PGF within the Dual Fermion approach. We 
will show that already a simple ladder approximation 
gives very reasonable results for the Hubbard model at 
half filling, for example a reduction of the critical tem- 
perature of the antiferromagnetic instability compared 
to the DMFT result is obtained already in the frame- 
work of the single-site calculations. It occurs at the Neel 
temperature which is in agreement with quantum Monte 
Carlo results 23 for the same parameters. It is remark- 
able that to obtain these results it is sufficient to take 
the bare dual fermion vertex as the irreducible vertex in 
the Bethe-Salpeter equations and only the lowest-order 
diagram for the dual self-energy. Such effectiveness is 
achieved by transforming the original interacting prob- 
lem to so-called dual fermion variables. By this we are 
able to include the local contribution to the self-energy 
into a bare propagator of the dual fermions and achieve a 
much faster convergence of the perturbation series. The 
outcome of the scheme is the 2PGF in the original vari- 
ables restored from the dual 2PGF with the help of an 
exact relation-^.. 

This paper is organized as follows: In Section II we 
review the dual fermion formalism in the general multi- 
band formulation. Section III gives a general overview of 
the exact equations for the 2PGFs in different channels. 
We then discuss various possible approximations which 
should be used for different purposes. We further derive 
the exact relation between 2PGFs in the dual and original 
variables. Section IV is devoted to the application of the 
approach to the two-dimensional (2D) Hubbard model 
at half filling and to the discussion of those results. In 
section V we give the conclusions and summary. 



II. DUAL FERMION FORMALISM 

The goal is to find an (approximate) solution to the 
general multiband problem described by the imaginary 
time action 

S[c*,c] 4k<rm + M) 1 - ^k ff ) mm / C wk(Tm / 

wkfrmm' 

+ J2 tfintte, Ci] . (1) 



Here h ka is the one-electron part of the Hamiltonian, 
u> n = (2n + l)%//3,n = 0, ±1, ... are the Matsubara fre- 
quencies, (3 and \i are the inverse temperature and chem- 
ical potential, respectively, a =f, J, labels the spin projec- 
tion, to, to' are orbital indices and c*,c are Grassmann 
variables. The index i labels the lattice sites and the 
k-vectors are quasimomenta. For applications it is im- 
portant to note that Hint can be any type of interac- 
tion. The only requirement is that it is local within the 
multiorbital atom or cluster. Consider, for example, the 



general Coulomb interaction 



p 

Hint [c* ,Ci} = 2 J dr [/1234c*! 



(2) 



where U is the general symmetrized Coulomb vertex and 
e.g. 1 = {wiTOiUi} comprehends frequency-, orbital- and 
spin degrees of freedom and summation over these states 
is implied. 

The idea of the dual fermion approach is to reformulate 
the lattice problem in terms of noninteracting impurities 
with their interaction replaced by a coupling to auxiliary, 
so-called dual fermions. 

As in the DMFT, we therefore introduce a local impu- 
rity problem in the form 

Simp[c*, c] = - C t<xm i( iuJ + M) 1 - A "<0mm' C ™' 

ujcrmm' 

+ H int [c*,c] , (3) 

where A is an as yet unspecified hybridization matrix 
describing the interaction of the impurity cluster with an 
electronic bath. 

Since we anticipate the decoupling of the interacting 
sites, we employ the locality of A to formally rewrite the 
original lattice problem in the following form: 

S[c , c] = ^ ] Sjmp [c^ j CTm , C u j am \ 
i 

C ojkerm ( A wc — h^o) mm i Cu>\i<rm> ■ 



(4) 



We introduce spinors c wk(T = (. . . , c^ kam , . . .), c* kcr = 
(. . . , c* kcrm , . . .). Omitting indices, the dual fermions 
are introduced via a Gaussian identity which in matrix- 
vector notation reads 

J exp (-f*Af -f*Bc-c*Bf)v[f*,f] = 



det(l) exp (c*bA~ 1 Be 



(5) 



This identity is valid for arbitrary complex matrices A 
and B. Choosing 



B = 9~l , 



(6) 



where g ua is the Green function matrix of the local im- 
purity problem in orbital space (to, to'), we obtain: 



5[c*,c,f*,f] = ^5 site ,i + 

i 

^ [Ck<7 9Zl (A^T - h ka ) 1 g~l fojker 



(7) 



Hence the coupling between sites is transferred to a cou- 
pling to the auxiliary fcrmions: 



Ssite,i — ^impfCj ? C^] 



^uicr 9uja c uia + C uia Q ua ^-u)ia ■ (8) 



Note that since g^a is local and the last term appears 
under a sum over all states labeled by k, the summation 
can be replaced by the equivalent summation over all 
sites. 

By transferring the inter-site coupling to interacting 
auxiliary fermions, only local degrees of freedom of the 
original fermions remain, so that we are able to integrate 
out the lattice fermions for each site separately: 



J exp (- S site [c* , Ci , £ * , fi] ) V [c| , a] 



(9) 



Formally this can be done up to all orders and in this 
sense the transformation to the dual fermions is exact. 
The above equation may be viewed as the defining equa- 
tion for the dual potential V[f*, f]. Expanding both sides 
of Eqn. ([9]) and equating the resulting expressions by 
order, one finds that the dual potential in lowest order 
approximation is given by 



V r [f*,f] = i^7S ) 34fafSfi4fi3 



(10) 



where 



7l234 — 9\\>9ti' 



imp 

Xl'2'3'4' 



imp,0 
Xl'2'3'4' 



ff3'3ff4'4 



imp.O 

Xl234 = 9l4ff23 - 313324 



(11) 



is the exact four-point reducible vertex for the original 
fermions which plays the role of the bare effective two- 
particle interaction for the dual fermions. The local two- 
particle Green function of the impurity model is defined 



as 



imp 

Xl234 — 



J irap 



ci c 2 Cg c\ exp ( - Simp [c* , c] ) V [c* , c] 



(12) 

The dual action now depends on dual variables only and 
can be written as 

s d [r,f] = -sc^Lr'w+E^.**] ■ ( 13 ) 

ukcr i 



The bare dual Green function is given by 



.9. 



(14) 



As we shall see below, for a properly chosen A, the 
DMFT result is already recovered within the simplest 
(zero-order) approximation for the dual potential V. In 
order to obtain the nonlocal corrections to the DMFT, we 
thus need to calculate the dual self-energy up to higher 
orders in V. This is achieved by performing a regular 




FIG. 1: The first two lowest order diagrams for the dual self 
energy E d . 



diagrammatic series expansion of the dual action, Eqn. 



We note that the only approximation is that in prac- 
tice the perturbation series expansion and the series for 
the dual potential need to be terminated at some point. 
Here we consider the first two lowest order skeleton di- 
agrams for Ed, constructed from the irreducible vertices 
and the dual Green function as lines. The use of skeleton 
diagrams ensures that the resulting theory is conserving 
according to the Baym-Kadanoff criterio n 20 ' 25 . The dia- 
grams considered here are shown in Fig. [TJ The lowest 
order diagram is local while the next diagram already 
gives a nonlocal contribution to the self energy. 

So far we have not established a condition for A, which 
was so far an arbitrary quantity. We require that the 
first diagram (Fig. QJi) in the expansion of the dual self- 
energy should be equal to zero at all frequencies. Since 
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(4) 



is local, we can use the condition ^ k = 0- In 
the simplest approximation, which corresponds to non- 
interacting dual fermions, the full dual Green function is 
replaced by the corresponding bare Green function and 
the above condition can be reduced to 



(A wcr — /ikcr) 



= 



(15) 



which is equivalent to the self-consistency condition for 
the hybridization function in DMFT. 

Thus, the case of non-interacting dual fermions corre- 
sponds to the full DMFT result for the original fermions. 
In order to go beyond DMFT it is sufficient to include 
a finite number of non- vanishing skeleton dual diagrams. 
In this paper all calculations are performed using only 
the lowest order non-vanishing skeleton diagram shown 
in Fig. QJ>). We postpone the explanation of the resulting 
numerical calculation procedure to section IV. 

The fact that we have employed an exact identity to 
transform to the dual variables has the important conse- 
quence that we can establish an exact relation between 
the lattice Green function and the dual Green function. 
This is important to obtain the Green functions for the 
original fcrmions without loss of information. We will 
further use this fact to establish the relation between 
original and dual 2PGFs. To this end, the partition func- 
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tion of the lattice is written in the two equivalent forms 

Z = J exp(-5[c*,c])D[c*,c] = 

Z f J J exp(-5[c*,c,f*,f])D[r,f]X>[c*,c] , 

(16) 

where 



Zf = JJ det [g ua (A wa - h^ a ) g ua \ 



(17) 



In the following, we introduce the quantities 

Luka = (A WCT - h^y 1 g~l (18) 

and 

Ruka = 9Zl - h^y 1 (19) 

which are matrices in orbital space. By taking the func- 
tional derivative of the partition function, Eqn. (|16[) . 
w.r.t. the Hamiltonian, the exact relation between the 
dual and lattice Green functions can now conveniently 
be written as 

Gujka — (A WCT — Tlker) + ^uku G^ k(T JJ^kcr (20) 

where the Green functions are defined via the imaginary 
time path integral as 

Gu = -~ /c 1 ^exp(-5[c*,c])X»[c* ) c] 

Zf J / 1 / 2 *exp(-5[f*,f])P[r,f] . 



G 



12 



z 



(21) 



III. TWO-PARTICLE GREEN FUNCTION. 

A. Relation between 2PGF's in the original and 
dual variables. 

Analogously to the single-particle Green function, the 
connection between the original and dual 2PGF's can be 
established by taking the second derivative with respect 
to h VfJ , of the partition function Z expressed in two differ- 
ent ways: Eqn. ([To]) . Here and further the small Greek 
letters denote the set {w,k, a, m} and the summation 
over repeated indices is implied. The 2PGFs are defined 
in a similar way as in Eqn. ()12[) : 

XX^ P = ^ J excite* exp(-S[c*,c])£>[c*,c] . (22) 



XU P = ^ J A/^/;«P(-5 d [f*,f])2?[f*,f] • (23) 



Varying the first expression for Z yields by definition 
1 6 2 Z 



Z Sh n \Sh 



— XXp,vp- 



(24) 



The most illustrative way to vary the second expression 
for Z from Eqn. (|16p is to write 



1 S 2 Z 1 S ( 5Z 



1 S 



ZSh 



Z 5h p \8h V n Z 6h p x \Sh u 



In the last expression we should use Eqn. ([20]) for G un- 
derstanding it as being diagonal in frequency, momenta 
and spin indices. From these two expressions we further 
obtain xx^ P = GxpG^ - ^f. Using the identity 



dA^ix) A -i dA iS A -i 



dx 



^ dx W 



for the derivative of an inverse matrix with respect to a 
parameter and the Eqn. (|20p we can rewrite the func- 
tional derivative of G as 

(A - [L G t RL„ + [L G„ fi] A „ (A - ft)" 1 .(25) 



The derivative in the last term in Eqn. (|25p is evalu- 
ated using the definition of the dual Green function, Eqn. 
(l2"Tj) and the identity 



Z 

Sh~x [Z~ f 



-Lxx 1 -^-G'kt n i R i, 



T X'p' ll p'p ■ 



(26) 



Thus, 

fry,* 
Sh p x 



Lxx' {Gx'p'Gp'v' - X\'n> v 'p>) R, 



p'p 



(27) 



Using Eqn. (|27|) in (1251) . we obtain the final result: 

XXfiup 



(A - hy 1 ® (a - hy 1 



(A-hy 1 ® [LG d R] 
[L Gd R]® {A-hy 1 



\p,vp 
\pvp 
\ P up 



+-C'AA / L^i Xx'ii'u'p 1 Ru'u Rp'p i (28) 

where (A ® B)x fiVp = AxpB^ - Ax v B lip is the anti- 
symmetrized direct matrix product. 

The formula ([28]) can be (symbolically) rewritten as 



X - G ® G = L L ( x d - G d <g> G d ) RR, 



(29) 



which shows quite clearly that the two particle excita- 
tions in the original and dual variables are identical. 
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B. Bethe-Salpeter Equations for Dual Fermions 

As is clearly seen from the end of the previous section, 
the problem of calculating the 2PGF in the original vari- 
ables reduces to the problem of calculating it for the dual 
fermions. In this section we review the standard Bethe- 
Salpeter equation approach to this problem and justify 
the approximations we use in the calculations, the details 
of the latter will be given in Section IV. 

In what follows we denote the 2PGF as x an d the ver- 
tices as T. Different sub- and superscripts will be used to 
distinguish between different channels and types of ver- 
tices. Small Greek letters stand for the set uj,k,m (or 
with appropriate changes u — > r and/or k — > x if the cal- 
culations are done in imaginary time or real space), spin 
will be taken into account explicitly. We also omit the 
index "d" for "dual" unless it is necessary for the sake of 
clarity. Otherwise all quantities used are dual by default. 

The full vertex V is defined by: 

crcr r~i<j s~ia j^aa f~ia s~iu 

XXpvp ^XX'^pfi' X' p.' v' p'^* v' v^- 1 p' p 

+GxpGp\,L> ~ Gx„G° p 5 aa i . (30) 

For this vertex one can write three different Bethe- 
Salpeter equations: 

pcrcr' _ -ppp,aa',ii , p y,pp,aa' ,in pa pa' perer' 
Xpvp 1 Xpup '^ Xpp'X' ^A'p' p'v'vpi 

(31) 

P<t<t ■pphl,crcr,ir -pphl, a, a, ir ^,cr s-ia -pcrcr 

1 Xpvp ~ 1 \fivp ~~ 1 Xi/'vX' ^A'p' p'f pVp'P' 

(32) 

yaa' -pphO,aa',ir _ y,phQ,aa",\r pa" pa" -pa" cr' 

Xpvp Xpvp Xp'X'p "AV^jiV "'w'' 

(33) 

where a denotes —a and £ := 1 — \& a a' ■ These equa- 
tions are written in different channels. The first one in 
the particle-particle channel and the latter two in the 
particle-hole channel. The two ph-channels differ from 
each other by the total spin (0,1) of the scattered particle- 
hole pair. r ir stands for the irreducible vertex in the 
corresponding channel. 

Obviously, the problem of calculating T is numerically 
solvable if r ir is known. The simplest possible approxi- 
mation for any irreducible vertex is the bare four-point 
interaction of dual fermions i.e. r )^\ In this paper we 
will not go beyond this approximation for practical cal- 
culations. Substituting it into Eqns. (|32l33p we obtain 
approximations for the exact vertex T in the respective 
channel: 

-ppp,aa' _ (4)oV * (4)oV pa pa' j^pp,aa' 
Xpvp 'Xpvp ' s IXpp'X' ^ X' p'^ p' v' L p'v'vpi 

(34) 

Xpvp ^Xpvp ^Xv' v\' \' p p' v p'pfi'p ' 

(35) 

-pphO.o-cr' (4)<t<t (4)(Tcr /-ht" /~icr" -pphO,a" a' 

Xpvp ^Xpvp ikp'X' p X' V*- 1 p' p' v'pvp' 

(36) 



We expect this approximation to capture the essential 
features of the problem if the channel is chosen prop- 
erly. Indeed, as can be seen from the previous works on 
dual fermions, the results for the Green function turn 
out to be qualitatively correct already in the lowest or- 
der in 7, this means that 7 is a good small parameter in 
the problem. Our expectations are justified by numerical 
calculations for the 2D Hubbard model at half-filling in 
the next section. It is worth mentioning here that we do 
not expect this simple approximation to give quantita- 
tively correct results for superconductivity, as the spin- 
spin correlations are known to be of high importance for 
the d-wave superconductivity of the Hubbard model. In 
order to approach this problem first the vertex in the 
ph-channel should be calculated and the result should be 
used as an approximation for the irreducible vertex in 
the pp-channel. This work is in progress now. 

We also want to comment on the standard DMFT ap- 
proach for calculating the 2PGFs. As has been men- 
tioned before, the DMFT result for Green function is 
reproduced from the dual fermion approach if the bare 
dual Green functions are used. On the other hand, if 
such an attempt is made for the 2PGF, i.e. the bare 
dual Green functions and the bare dual four-point vertex 
7 are inserted in the Eqn. (|28[> . we do not reproduce the 
DMFT result. Instead what we obtain is: 

X-G D ®G D = LL(G d ' G d '°^G d ' G d >°)RR , (37) 

where G D = (g- 1 + (A-h))~ 1 is the DMFT lattice Green 
function. Using the relation (A - h)~ 1 g- 1 G d '° = -G D 
(which can be proven by straightforward calculation) 
leads to: 

X-G D ®G D = G D G D 7 ^G D G D (38) 

In the last expression 7W - being the full vertex of the 
Anderson impurity model - can be symbolically written 
in form of a Bethe-Salpeter equation: 

7W = 7 ir +7 ir fffl7 (4)_ (39) 

If here we replace the local Green functions g in the lad- 
der by G D , we immediately restore the conventional re- 
sult for the 2PGF in DMFT approach: x = Xo + Xol"x 
with xo = G D <E> G D . But this substitution can only 
be done "by hand" and as is obvious from the structure 
of the dual diagrams cannot be reproduced automati- 
cally at any order of the perturbation theory. This is 
an indication that the conventional way to calculate the 
2PGF is not <I>-derivable and does not constitute a sys- 
tematic way to calculate the fc-dependent vertex function. 
This drawback is excluded in our approach. However, 
the DMFT approach still works reasonably. One can see 
that by noticing that g = G D — G d '° (derived straightfor- 
wardly). The local part of G d, ° vanishes by construction 
and the non-local part is generally small. Thus one can 
say that the DMFT approach to calculating the 2PGF 
corresponds to considering the first non-trivial diagram 
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for the dual 2PGF with a small change of a part of the 
internal Green functions of the vertex. 

Finally we would like to mention that apart from 7^ 
also higher cumulants are present in the dual fermion 
approach. We do not expect them to be important for 
calculating the 2PGF and do not take them into account 
in this paper. However, the 6-point vertex 7W which 
plays the role of the effective 3-particle interaction for the 
dual variables can turn out to be of crucial importance 
for non- linear optics calculations. 



IV. APPLICATION TO THE 2D HUBBARD 
MODEL 

Let us now turn to the Hubbard model. The 2D Hub- 
bard model is described by the Hamiltonian 



H 



E 



(40) 



A ^ > 9,1 W 



DMFT 



G d 



outer loop 



inner 
loop J, 



A new = A old+5 - 1 Gf oc G 1 -i 



FIG. 2: (Color online) Illustration of the dual fermion calcu- 
lation procedure. 



with the bare dispersion relation given by = 
—2t(cosk x + cosfcy). In the following the energy scale is 
fixed by setting the nearest neighbor hopping iy = t = 1. 
In this study, we restrict ourselves to the case of half- 
filling and postpone the case of finite doping including 
an analysis of the superconducting instability to a later 
publication. 

It is known that at low temperatures strong anti- 
ferromagnetic correlations develop in the 2D Hubbard 
model at half-filling. This is due to the perfect nest- 
ing of the Fermi surface. Rigorously speaking, a tran- 
sition to a long-range ordered state at finite tempera- 
tures cannot occur since this would break the continu- 
ous spin symmetry, which is prohibited by the Mermin- 
Wagner theorem^. For calculations however, one should 
keep in mind that a transition to the antiferromagnetic 
state is possible: Being bound to approximations, the 
fluctuations responsible for destroying the long-range or- 
der are not fully accounted for and the implications of 
the Mermin- Wagner theorem are not applicable. In this 
section we study how the antiferromagnetic instability 
emerges within our approach and how the transition tem- 
perature compares with that obtained within the DMFT. 

The calculations are performed for the 2D Hubbard 
model at half-filling within the paramagnetic phase. In 
this publication, we restrict ourselves to the application 
of the single-site dual fermion approach. As a prereq- 
uisite for a better understanding on how the corrections 
to the DMFT emerge, we first briefly review the calcula- 
tion procedure. Then, in order to illustrate the method 
and to underline some of the implications from the sus- 
ceptibility calculations, we first present some results for 
single-particle properties obtained from DMFT and dual 
fermion (DF) calculations. Afterwards, explicit expres- 
sions for the solution of the Bethe-Salpeter equation and 
the calculation of susceptibilities are introduced. Then 
we discuss the two-particle properties obtained using the 



ladder approximations to the 2PGF discussed in the pre- 
ceding section. In order to compare the DF and DMFT 
results, we transform the dual susceptibility to the corre- 
sponding result for the original lattice fermions. The DF 
results are then contrasted to the corresponding DMFT 
results for different values of the on-site repulsion and 
different temperatures. We find considerable corrections 
to the DMFT from our DF calculations, in particular fol- 
iar ge values of U. 



A. Dual Fermion calculations 

Since the method is rather new, let us first briefly re- 
view the calculation procedure. Additional information 
can be found in Refs. HO and Hi]. Each calculation is 
started with a regular self-consistent DMFT calculation. 
This is achieved by requiring the local part of the bare 
dual fermion propagator (i.e. no corrections to the dual 
self-energy are taken into account as indicated by the 
crossed diagram in Fig. [2]) to be zero. We then calculate 
the vertex . This enables us to calculate an approx- 
imation to the dual self-energy by summing up the first 
diagrams in the perturbation series expansion. From this 
and the bare dual Green function an approximation to 
the dual Green function G d is obtained via Dyson's equa- 
tion, which is subsequently used in the diagrams. This 
inner loop is executed until self-consistency. This effec- 
tively replaces the diagrams depicted in FigfT]by the cor- 
responding skeleton diagrams. In general, the first-order 
diagram is now non-zero, in contrast to the DMFT result. 
From the condition that this diagram should be zero, 
we construct a new hybridization function, which serves 
as input for the calculation of a new local Green func- 
tion and renormalized vertex 7^ in the impurity solver 
step. We repeat this outer loop until self-consistency 
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is reached. The computational cost for the calculations 
aside from DMFT is comparable to the DMFT itself, 
whereby the calculation of the vertex is the computation- 
ally most expensive part. We use the continuous-time 
quantum Monte Carlo impurity solve r 27 i 28 i 29 for the so- 
lution of the impurity problem and for the calculation of 
the vertex. 



B. Single- particle properties 

In this section we present some characteristic single- 
particle properties to illustrate our approach. All calcu- 
lations have been performed for U = 4 and 8. Note that 
the bandwidth is W = 8t = 8. 

To begin with, we present results for the local density 
of states (see Fig. [3]). For U — 4, we find qualitatively 
the same behavior within the DMFT and DF calcula- 
tions, with the difference in the DF result being a sup- 
pression of the spectral weight around the Fermi energy 
compared to DMFT. For this value of U the gap does not 
open, but in should be noted that in order to fully open 
a gap e.g. within the dynamical cluster approximation 
(DCA) requires large clusters of the order of 64 site s 30 ' 31 . 
At U = 8 however, we find a drastic change in the lo- 
cal density of states when the dual corrections are taken 
into account. We observe qualitatively similar changes 
in a broad temperature range above the critical temper- 
ature. In the DF calculations, the spectral weight around 
the Fermi level is strongly suppressed. This results in a 
pseudogap which persists up to higher temperatures. It 
is believed to be caused by non-local spin fluctuations, 
which are not present in the DMFT. 

This further suggests that non-local corrections to the 
DMFT self-energy become important in this parameter 
regime. To underline this picture, we show the real and 
imaginary parts of the k-dependent lattice self-energy on 
the first Matsubara frequency. Within DMFT, this quan- 
tity is just a constant, £(7rT, k) = S(7rT). For the case of 
large on-site repulsion, we find a strong renormalization 
of the bare dispersion law hu as shown in Fig. 0] Also 
the imaginary part of the self-energy exhibits a strong 
k-dependence (see Fig |S) . Qualitatively similar features 
are found for higher temperatures. 

A quantity which is directly accessible to experimen- 
tal observation via angular resolved photo emission spec- 





FIG. 3: Local density of states for the 2D Hubbard model 
at half-filling obtained within DMFT (dashed lines) and dual 
fermion calculations (solid lines) for U/t = 4 (left) and U/t = 
8 (right) at inverse temperature j3 = 4.5. 




FIG. 4: (Color online) Contour plot of the real part of the 
lattice self-energy on the first Matsubara frequency 5ft£(u; = 
vrT, k) (centered at the T-point) for U/t = 4 (left) and U/t = 8 
(right), calculated on a grid of 256 x 256 k-points at inverse 
temperature f3 = 3.5. 
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FIG. 5: (Color online) Contour plot of the imaginary part 
of the lattice self-energy 9E(lj = nT, k) (centered at the F- 
point) for U/t = 4 (left) and U/t — 8 (right) at inverse tem- 
perature (5 — 3.5. 



troscopy is the single-particle spectral function A(k,oj). 
In Fig. [6] we show A(k, w) for {7 = 8 and /3 = 4.5 ob- 
tained from the DMFT and DF calculations, respectively. 
As expected from the local density of states, the DMFT 
and DF spectral functions are quite different. The most 
striking feature is the absence of a coherent peak at the 
Fermi level in the DF calculations. Our DF spectral func- 
tion very well resembles the characteristic features of the 
one given in Ref. In this work, an approach based on 
Hubbard operators including spin-fluctuation corrections 
to the self-energy has been used. The self-energy was 





FIG. 6: Spectral function for the 2D Hubbard model at half- 
filling obtained within DMFT (left) and DF (right) calcula- 
tions for U/t — 8 at inverse temperature /3 = 4.5. From bot- 
tom to top, the curves are plotted along the high-symmetry 
lines r —* X — ► M — * T. The high symmetry points X and 
M are marked by dashed lines. 
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FIG. 7: (Color online) Comparison of the DMFT Weiss field 
go(ioj) — [iw + fi — A(iaj)] -1 and its modification after a fully 
self-consistent dual fermion calculation (DF), for U/t = 4 
(left) and U/t — 8 (right). The dashed curves correspond 
to the inverse temperature /3 = 3.5 and the solid curves to 
f3 = 4.5. 



obtained from the self-consistent solution of a two-site 
dynamical cluster problem and is thus computationally 
more involved than our single-site approach. In addition, 
it treats the local (i.e. nearest-neighbor) correlations ex- 
plicitly. This option is also present within the cluster- 
formulation of our approach. 

Since obviously the DMFT insufficiently describes the 
physics leading to the appearance of the pseudogap, one 
might expect that the dynamical mean field, which de- 
scribes the influence of the surrounding sites in a mean 
field manner, does not accurately reflect the surrounding 
in the real system. In such a case, the renormalization 
of the hybridization function A(icj) (cf. Fig. [2]) should 
become important. This can be seen in Fig. where we 
compare the Weiss field go(iu>) — [iui + fi — A(iw)]~ 1 from 
the DMFT and DF calculations. For U = 4, the change 
is small and almost does not depend on temperature. On 
the other hand, for U — 8, we find a strong renormaliza- 
tion, which significantly increases when the temperature 
is lowered. This reflects the fact that the local singlet for- 
mation prevails at low temperatures. Also note that the 
DMFT result hardly depends on temperature, in contrast 
to the DF results. 



C. Solution of the Bethe-Salpeter equation 

In order to calculate the susceptibilities, we need to 
calculate the reducible two-particle vertex by solving the 
Bethe-Salpeter equation. Explicitly, the Bethe-Salpeter 
equations in the singlet and triplet particle-hole channels 
(see Fig. [8]) are 

-pph0,aa' i \ _ oV ^ \ * \ * -,aa" s-i&cr" 

xCi(k + q)rSV'(q) 

(41) 



P u" k 

xC + o(k + q)r^(q). 

(42) 

Here 7 is the irreducible vertex and G^ CT (k) denotes the 
fully self-consistent nonlocal dual fermion propagator. To 
be complete, we wrote the spin indices on the propaga- 
tors. The calculations however are carried out for the 
paramagnetic case, i.e. G da = G dcr . Capital letters 
denote the bosonic and small letters fermionic Matsub- 
ara frequencies. N is the linear dimension of the lattice 
and d the dimension. Due to the ladder approximation, 
the resulting fully reducible vertex T only depends on a 
single momentum q. 

The Bethe-Salpeter equation is solved iteratively. A 
starting guess for the fully reducible vertex T — is in- 
serted into the equation and a new approximation rw is 
obtained. This process is repeated until self-consistency. 
As a convergence criterion we require the difference be- 
tween two successive iterations to be smaller than some 
predetermined accuracy e, i.e ||r( n+1 ) — rw|| < e. As 
a measure for the deviation we use the entrywise norm 
||A|| = J2ij\ a ij\- I n principle it is possible to solve the 
Bethe-Salpeter equation by supermatrix inversion (which 
requires a decoupling of the singlet channel into spin and 
charge channels). However, close to the instability one 
would be faced with the intricate task of inverting ill- 
conditioned matrices. When iterating the equation the 
instability is signaled by a decelerated convergence. This 
is related to the fact that the leading eigenvalue of the 
corresponding matrix tends to one, so that more and 
more diagrams in the ladder need to be taken into ac- 
count in order to obtain an accurate result. If one is 
not interested in the susceptibility itself but only in the 
instability, this can be circumvented by locating the in- 
stability by finding the parameters for which the lead- 
ing eigenvalue of the BSE-related eigenvalue problem be- 
comes one. 

Note that the convolution of the two Green functions 
in Eqns. (|41I42[) involves a sum over n = N d k- vectors for 
each value of q and thus of the order of 0(n 2 ) arithmetic 
operations. This becomes tedious already for relatively 
small lattice sizes. It is then possible to reduce the order 
down to O(nlogn) by performing a so-called fast convo- 
lution using fast Fourier transforms. The convolution is 
calculated for all vectors q simultaneously and stored in 
memory. It is also used for the calculation of the suscep- 
tibilities (see below). Efficient transforms are provided 
by standard packages^. 

Once we have found a converged solution for the two- 
particle vertex T, we are able to calculate nonlocal sus- 
ceptibilities. For the paramagnetic case, we have 

(5,^)(n > q3 = i«n T n T )-(n T n i ))(n j q) , (43) 

where (n o .n CT /)(0,q) = Xo ff '(fyq) + X aa ' i^, q)- This is 
valid for lattice and dual fermions. For the dual fermions, 
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FIG. 8: The Bethe Salpeter equation in the three different 
channels. In the calculations we considered the equations for 
the particle-hole channel. 



the first term is given by the bubble diagram (dual sus- 
ceptibilities are marked by the tilde): 

^'(^i) = -^EE G ' ff ( k )C( k + q) , (44) 



uj k 



while the nontrivial part of the susceptibility, cf. Fig. [9l 
is given by 

r <T '(^,q) = 



f3 2 N 2d 



kk' 



xGj?'(k')G^ n (k' + q) 



(45) 

For the other channel, the nontrivial part of the suscep- 
tibilities (S ± S^) is given by {c\c- a c\c a ) {ft, q) = X q 9 + 
X aa ■ For the dual fermions,this is obtained via 

X CTff (tt,q) = 
1 



p2 N 2d 



EEC'(q)C(k)C(k + q) >< 

xG^(k')G^ +f2 (k'+q) . 



uju)' kk' 



(46) 

In order to be able to compare with known results, 
we need to transform the dual susceptibility to the sus- 
ceptibility for the original fermions. According to Eqn. 
(p?9")) . this can be achieved by first multiplying the ver- 
tex by four dual Green function legs to obtain the dual 
two-particle Green function and then multiplying by the 
functions L and R defined in Eqns. (|18I19|) . In order 
to calculate the susceptibility (S z S Z )(Q, q) this merely 
amounts to calculate the susceptibility according to Eqns. 
(I43H45|) , but with the dual Green functions in Eqn. (|4"5|) 
replaced by the modified propagators 

Q L = (A WCT - /ikcr) -1 9Za G ta 



G ta 9ujl 



(47) 
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FIG. 9: Diagrams for the susceptibilities in the two different 
particle-hole channels 



(these are identical for a diagonal basis) and the Green 
functions in Eqn. (|4~4"|) by the lattice Green function. 

The DMFT susceptibility can be straightforwardly ob- 
tained since we already have the fully reducible vertex of 
the local impurity problem, 7", at our disposal. From 
this we can obtain the corresponding vertex in the spin 
channel as = 7+t — 7+f • The irreducible vertex of 
the impurity problem is then obtained by matrix inver- 
sion using the standard relation 



(48) 

= £ k G5(k) is 



where xln = ~^^G^ and G^ loc - 
the local DMFT Green function. The vertex is obtained 
by iterating the BSE 



)=7i' 



1 



(j\rd z2 /2 7sww"n G °' ( k ) x 

u>" k 

xG u D „ +S2 (k+q)r» V o(q) . 



(49) 

where G° is the DMFT lattice Green function, as before. 
The susceptibility itself is obtained using equations simi- 
lar to Eqns. (|44I46|) with the vertex and Green functions 
replaced by the appropriate DMFT quantities. 



D. Two- particle properties 

Now we turn to the two-particle properties. In order to 
illustrate the statement that two-particle excitations are 
the same for real and dual fermions, we first consider the 
dual susceptibility, Eqn. PS)) . In Fig. [TUl we show the 
results for the nontrivial part \ zz °f the dual susceptibil- 
ity as a function of temperature for U = 4 and different 
k-points. This calculation was performed for an 8 x 8 lat- 
tice. The susceptibility diverges at the wave vector (w, 7r), 
indicating a transition to the antifcrromagnetic ordered 
state, while at the other k-points it does not show any in- 
dication for a divergence. The susceptibility of the lattice 
fermions diverges at the same temperature, as expected. 
We would like to note that the dual bubble diagram, Eqn. 
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FIG. 10: (Color online) Nontrivial part of the dual suscepti- 
bility \ zz = — X TJ ") f° r different k-points as a function 
of temperature. The divergence at q = (n, n) indicates the 
antiferromagnetic instability. 



FIG. 11: (Color online) The antiferromagnetic correlations 
visualized in real space for = 4.25 and f3 = 4.5. At temper- 
atures close to the transition, the correlations become essen- 
tially independent of distance. 



(|44| can be negative so that the dual susceptibility Xo + X 
becomes negative at high temperatures. This is due to 
the fact that the spectrum of the dual Green function 
need not be positive-semidefinite, which follows from the 
fact that its local part exactly vanishes. Indeed, since the 
dual fermions represent the nonlocal part of the lattice 
fcrmions, they are no physical particles. However, this 
does not affect the lattice susceptibility as shown below 
and we did not encounter any non-analyticity problems 
in all our calculations. 

In the paramagnetic case the relation ^((S + S~) + 
(S~S + )) = 2{S Z S Z ) holds. We find that this relation is 
very well fulfilled within our calculations. It follows from 
the invariance under rotations in spin space and it can 
be generally proven 34 that r~ := T*™ a<T - T P J% S (which 
describes the interactions of the particle-hole pair in the 
triplet state with spin projection o z = 0) is equal to T phl 
(which corresponds to the spin projections a z — ±1, de- 
pending on the sign of a). This identity is readily shown 
to be preserved by the BSEs Eqns. (|41I42[) given that 

IcrSScr = laaaa ~ laaaa holds, which is fulfilled in Our 




FIG. 12: (Color online) Momentum dependence of the static 
susceptibility (centered at the M-point) for (7 = 4 and differ- 
ent temperatures calculated on a 256 x 256 lattice. Top row 
from left to right: f3 = 0.5, (3 = 1.0 and /3 = 3.0. Bottom row: 
P = 4.0, P = 4.5 and (3 = 4.65. 



calculations up to a small numerical error. 

Using the definition of 7^ 4 \ Eqn. (|TT|l ■ one finds that 
this small numerical difference can be traced back to the 
error due to Monte-Carlo (MC) averaging of the two- 
particle Green function: In each MC measurement, one 
measures a quantity G12 which corresponds to the con- 
tribution of the particular configuration to Green's func- 
tion (for details, see Ref. I^glll). The Green func- 
tion G12 = (G12) is obtained as the MC average of this 
quantity. Similarly, using Wick's theorem, the two par- 
ticle Green function X1234' 7 i s obtained as (G^G?^) — 



Thus one has XiSu 



aaaa 
A1234 



{G\ A G%2 

Xl234 — 1^12^34 



-<G? 4 Gf 2 } and 



(G14G32)) - (Gi 2 G^ 4 ) 



Since the quantities G12 for different spins can differ even 
for paramagnetic systems, the above quantities can only 
be equal within the MC error. 

In Fig. QTJ we show the real space dependence of the 
spin-spin susceptibility for the 8x8 lattice. At lower 
temperatures, the correlations decay quickly as a func- 
tion of distance, although the antiferromagnetic charac- 
ter is clearly visible. At lower temperatures close to the 
transition, the correlation length approaches the lattice 
size, and the susceptibility becomes essentially constant 
as a function of distance. 

In Fig. [T^lwe illustrate the k-dependence of the suscep- 
tibility (S Z S Z )(Q = 0, k). For low temperatures, the sus- 
ceptibility appears to be delocalized in k-space, leading 
to the fast decaying correlations as a function of distance. 
However, even for the highest temperature the maximum 
is located at k = (tv,tt), exposing the tendency to anti- 
ferromagnetism. For low temperatures, the susceptibility 
becomes strongly peaked in the (n, tt) direction. 

As far as the single-particle properties are concerned, 
our approach has proven to give physically correct re- 
sults. However, since the calculation of the susceptibil- 
ities relies on an additional approximation, namely that 
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FIG. 13: (Color online) Inverse of the antiferromagnetic sus- 
ceptibility xaf = (S z S z )(fl = 0, q = (7r,7r)) as a function 
of temperature as obtained from DMFT and within the dual 
fermion approach for (7 = 4 and different lattice sizes. Here 
e.g. the label 64 x 64 indicates the number of k-points used in 
the Brillouin zone integration. The inset compares the fully 
self-consistent DF result (dashed line) with the one obtained 
from the first iteration of the outer loop (solid line) shown in 
Fig.H 
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FIG. 14: (Color online) The antiferromagnetic susceptibility 
for U — 8. The inset shows the same on a larger scale. For 
high temperatures, the dual fermion result converges to the 
DMFT result. 



we take the bare interaction of the dual fermions, 7^ 4 - ) , as 
an approximation for the irreducible vertex, we need to 
demonstrate that the approach still yields sensible results 
for the two-particle properties. 

A stringent test for this approximation is whether the 
results are still an improvement compared to the DMFT. 
In Figs. 1131141 we therefore plot the inverse antiferro- 
magnetic susceptibility xaf = (S Z S Z )(^ = 0,k= (7r,7r)) 
as a function of temperature. The critical temperature 
is given by the point to which the corresponding curve 



extrapolates to zero. Note that the exact solution is ex- 
pected to have a Neel temperature Tjv = 0. One might 
expect that since the DF and DMFT results are rather 
similar for U = 4, this also applies for the susceptibil- 
ity. Still, we find that the DMFT result is systematically 
below the DF result, and the DF critical temperature is 
correspondingly lower than in the DMFT. This can be 
interpreted as an effect of the non-local spin fluctuations 
which are not accounted for in the DMFT and which sup- 
press the critical temperature. In order to understand 
the scaling behavior one should recognize that unlike the 
DCA, which essentially is an expansion in 1/N C (N c is 
the clustersize) and hence becomes formally exact in the 
limit of infinite cluster o 35 ' 36 , our approach does not be- 
come exact in the limit of infinite lattice sizes. Thus 
we cannot reproduce the scaling behavior as observed in 
DCA calculations 37 . This is due to the termination of the 
perturbation series. On the other hand, our approach be- 
comes formally exact if all diagrams are considered in the 
limit of infinite lattice sizes. Since we sum a perturbation 
series for auxiliary fermions, an a priori statement which 
diagrams include what kind of physics in terms of the lat- 
tice fermions is not possible. However, we find that the 
approach is essentially converged for lattices larger than 
16 x 16, setting a scale for the range of the fluctuations 
included. We expect this scaling to be more significant 
if either the number of diagrams considered is increased, 
or the starting point is improved, i.e. the dual fermion 
calculation is performed on top of an e.g. 2x2 cluster- 
DMFT calculation. 

The inset of Fig. Q2] compares the antiferromagnetic 
susceptibility obtained from a fully self-consistent DF cal- 
culation (dashed line, with diamonds) to that obtained 
from the first iteration of the outer loop, i.e. just one 
inner loop as shown in Fig. [2] has been performed. The 
hybridization is thus the one obtained from DMFT. One 
sees that the fully converged DF result has a slightly 
larger critical temperature than the intermediate result. 
We have observed this tendency in all our calculations. 
We attribute this to the renormalization of the hybridiza- 
tion as compared to DMFT, which pronounces the effect 
of the local singlet formation and thus favors antiferro- 
magnetism. 

In Fig. Q3] we show the temperature scaling of the in- 
verse susceptibility for U — 8. Here we find that almost 
no scaling with the lattice size is visible. This can be 
attributed to the fact that the physics of singlet forma- 
tion prevails in this parameter regime and this is what 
is primarily mediated by the first diagram in the pertur- 
bation expansion. This is consistent with the fact that 
the transition temperature is higher compared to U = 4. 
We thus expect a considerable reduction of the transition 
temperature when the local singlet-formation is treated 
explicitly within the cluster-extension of our approach. 
However, the fluctuations are still taken into account and 
lead to a suppression of the critical temperature which is 
somewhat more pronounced than for [7 = 4. 
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V. CONCLUSIONS 

To conclude, we have proposed a scheme to calculate 
the two-particle Green function within the dual fermion 
framework. This enables us to study two-particle interac- 
tions in strongly correlated systems where non-local cor- 
relations cannot be neglected. We formulated the Bcthc- 
Salpeter equations for the full vertex in the particle- 
particle and particle-hole channels and proposed an ap- 
proximation for practical calculations. A possible exten- 
sion of this scheme has been proposed which allows the 
application to superconductivity while taking also the an- 
tiferromagnetic fluctuations into account. We have estab- 
lished an exact relation between the two-particle Green 
functions in dual and original variables, which ensures 
that two-particle excitations of real and dual variables 
are identical. The identity was used to transform the 
dual susceptibility to the one of the original fermions. 
Within our proposed approximation we have applied the 
scheme to the 2D Hubbard model at half filling in the 
paramagnetic phase. We find strong modifications of 
the single-particle properties in our dual fermion calcula- 



tions compared to the DMFT. This is encouraging given 
that our calculations were performed within the single- 
site formalism. Concerning the two-particle properties, 
we also found an improvement compared to the DMFT. 
The critical temperature is suppressed due to the incor- 
poration of the non-local spin fluctuations, thus show- 
ing that this approximation goes well beyond the con- 
ventional DMFT scheme for calculating the two-particle 
Green function. We expect a further improvement of the 
results by the cluster formulation of our approach, which 
explicitly takes the local singlet formation into account. 
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